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ABSTRACT 

We present a new, fast, algorithm for the separation of astrophysical components su- 
perposed in maps of the sky. The algorithm, based on the Independent Component 
Analysis (ICA) technique, is aimed at recovering both the spatial pattern and the 
frequency scalings of the emissions from statistically independent astrophysical pro- 
cesses, present along the line-of-sight, from multi-frequency observations, without any 
a priori assumption on properties of the components to be separated, except that all 
of them, but at most one, must have non-Gaussian distributions. 

The analysis starts from very simple toy-models of the sky emission in order to 
assess the quality of the reconstruction when inputs are well known and controlled. In 
particular we study the dependence of the results of separation conducted on and off 
the Galactic plane independently, showing that optimal separation is achieved for sky 
regions where components are smoothly distributed. 

Then we move to more realistic applications on simulated observations of the 
microwave sky with angular resolution and instrumental noise at the mean nominal 
levels for the Planck satellite. We consider several Planck observation channels con- 
taining the most important known diffuse signals: the Cosmic Microwave Background 
(CMB), Galactic synchrotron, dust and free- free emissions. A method for calibrating 
the reconstructed maps of each component at each frequency has been devised. The 
spatial pattern of all the components have been recovered on all scales probed by the 
instrument. In particular, the CMB angular power spectra is recovered at the percent 
level up to Imax — 2000. 

Frequency scalings and normalization have been recovered with better than 1% 
precision for all the components at frequencies and in sky regions where their signal- 
to-noise ratio ^ 1.5; the error increases at ~ 10% level for signal-to-noise ratios ~ 1. 

Runs have been performed on a Pentium III 600 MHz computer; although the 
computing time slightly depends on the number of channels and components to be 
recovered, FastICA typically took about 10 minutes for all-sky simulations with 3.5' 
pixel size. 

Although the quoted results have been obtained under a number of simplifying 
assumptions, we conclude that FastICA is an extremely promising technique for 
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analyzing the maps that will be obtained by the forthcoming high resolution CMB 
experiments. 

Key words: methods - data analysis - techniques: image processing - cosmic mi- 
crowave background. 



1 INTRODUCTION 

NASA's Microwave Anisotropy Probe (MAP) satellite (Wright 1999) has begun its mission aimed at all-sky imaging of the 
last scattering surface of CMB photons, at several frequencies and with high angular resolution. In 2007, the very early phase 
of our Universe will be probed with even higher sensitivity and angular resolution by ESA's Planck satellite (Mandolesi et al. 
1998; Puget et al. 1998). In the meantime, new ground-based and balloon-borne experiments will provide improved maps of 
regions of the sky (see de Bernardis et al. 2001 and references therein). The main goal of these projects is to extract a map as 
accurate as possible of CMB anisotropies at all scales relevant for cosmology. This will allow us to determine the fundamental 
cosmological parameters with very high accuracy and to accurately test the models of structure formation. 

Maps produced by these missions will contain not only the CMB signal but also contributions from different components 
placed between us and the last scattering surface. These consist of Galactic free-free, thermal and spinning dust, and syn- 
chrotron emissions, of compact galactic and extragalactic sources as well as of the thermal and kinematic Sunyaev-Zeldovich 
effect in cluster of galaxies. Therefore, in order to exploit the wealth of cosmological and astrophysical information encoded 
in these maps, it is necessary to separate the different components with high accuracy and reliability. 

A great deal of work in this direction has been carried out by several authors (de Oliveira-Costa & Tegmark 1999, Tenorio 
et al. 1999, Hobson et al. 1998, 1999, Stolyarov et al. 2001, Prunet et al. 2001 and Baccigalupi et al. 2000) exploring different 
reconstruction techniques, from Wiener filtering and Maximum Entropy Method (MEM), both in real and Fourier space, to 
wavelets and neural networks. 

Traditional methods (Wiener filtering and MEM) have been successfully applied to all-sky maps. They need the prior 
knowledge of the frequency dependence of the astrophysical components to be recovered, as well as of the power spectra of their 
spatial pattern. Recently Stolyarov et al. (2001) have achieved a high quality reconstruction of the various components using 
a MEM algorithm working in harmonic space, assuming prior knowledge of the frequency dependence and of the azimuthally- 
averaged power spectrum of each input component except for the CMB. Computational requirements are quite demanding: 
the reconstruction requires 6 hours using 30 R10000 processors on an SGI Origin 2000 supercomputer and 14 Gb of RAM. 

The Independent Component Analysis (ICA), based on a feed- forward neural network (Baccigalupi et al. 2000), works 
remarkably well without any a-priori knowledge on the components, provided that they are statistically independent and have 
all, but at most one, non-Gaussian distributions, as is actually the case since, at most, the distribution of CMB fluctuations 
is Gaussian, while foregrounds are all highly non-Gaussian. The main advantage of this approach is that the algorithm is 
able to learn how to reconstruct independent components and their frequency scalings from the input maps. This technique 
is very interesting for at least three reasons. First, as we show in this work, it is able to recover independent components and 
their frequency dependencies with high accuracy also in the presence of instrumental noise. Second, it can be used as prior 
estimator for non-blind separation algorithms. Third, a blind approach is necessary whenever reliable priors for the foreground 
emissions are not available, as in the case of CMB polarization. Baccigalupi et al. (2000) have applied this technique to small 
sky patches of 15° x 15°, without instrumental noise. In these ideal conditions, CMB is reconstructed with 1% accuracy, 
weaker galactic components are recovered with ~ 10% error, while point sources are extracted down to a flux corresponding 
to ~ 0.7<TCMB- 

We present here the application of a new ICA method (FastICA), proposed by Hyvarinen et al. (1998), to the problem 
of separating the various astrophysical components present in all-sky maps. This technique proved to be extremely fast: on 
a Pentium III 600 MHz PC, it took only ~ 10 minutes to perform the component separation over the whole sky with ~ 3'5 
resolution. Furthermore, and importantly, this technique allows us to deal also with the instrumental noise, at least as long 
as it can be approximated by a Gaussian distribution and is uniformly distributed over the sky. 

This paper is organized as follows. In Section^ we formalize the problem of component separation and list the assumptions 
we have made. Section ^ describes in detail the FastICA technique used in this work. In Section ^ we study the quality 
of signal reconstruction with FastICA for different kinds of toy models i.e. when the input signals are well known and 
controlled. Section ^ deals with more realistic simulations of sky maps, focusing on the Planck satellite frequency channels 
and instrumental characteristics (noise and FWHM). Finally, a discussion and conclusions are presented in Section []. 
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2 FORMALIZATION OF THE SOURCE SEPARATION PROBLEM 

Let us consider the sky radiation as a function of direction f and frequency v, and assume that it results from the superposition 
of signals coming from jV different physical processes. 

N 

x{r,v)= y ^2,s j (r,v) (1) 

3 = 1 

where x is the total radiation from direction f at frequency v, and Sj are the individual physical processes that form the total 
radiation. 

The signal x is observed through a telescope, whose beam pattern can be modelled, at each frequency, as a shift-invariant 
point spread function B(f,v). For each value of v, the telescope de-focuses the physical radiation map by convolving it with 
the kernel B. The frequency-dependent convolved signal is input to an Af-channel measuring instrument, which integrates 
the signal over frequency on each of its channels and adds some noise to its output. The final signal at the output of the 
measurement channel at frequency v is given by 

f r 

3=1 

where t u (y') is the frequency response of the channel, and £ v {f) is the related noise map. To simplify the data model in 
Eq. (Q), let us now make the following assumptions: 
Assumption 1 

Each source signal is a separable function of location and frequency: 

«,-(<» = «i(f)/j(") j = l,...,N (3) 
Assumption 2 

B(r, v) is frequency-independent over each channel pass-band i.e. for each frequency where U(y) 7^ in the i-th channel, the 
corresponding beam function Bi is the same. These two assumptions lead us to a new data model: 

JV 



J2B v (r)*s 3 (f)- J t v {v')U{v')dv +e„(f) 
3=1 



= B v {?)*2^a VJ s 3 {r) + e v {r) (4) 
3=1 

where * denotes convolution, and 

a v] = J t v {u')fj{v')du' . (5) 

This data model can be further simplified by making another assumption: 
Assumption 3 

The radiation pattern of the telescope is the same for all the measurement channels, that is, the function B(f, v) is frequency 
independent: 

B(f,v) = B(f). (6) 
In this case, Eq. (^]) can be rewritten in vector form as: 

x(f) = As(r) * B(f) + e(f) = As(f) + e(f) (7) 

where each component, Sj, of the vector s is given by the corresponding source function Sj, as defined in Eq. (Eh, convolved 
with the frequency-independent radiation pattern B. The matrix A is the mixing matrix whose elements are the a v j defined 
by Eq. §. 

We detailed this derivation to stress all the simplifying assumptions made. Depending on the particular astrophysical 
application, one or more of these assumptions might be not justified. Each of these cases needs to be addressed by a specific 
strategy. In the present simple case, we are just assuming that the noise is additive, signal-independent, white, Gaussian, 
and stationary. For example, it can be noted that, in principle, assumption 3 is somewhat bad for applications where the 
instrument has different B v at different frequencies, since it requires that higher resolution channels are degraded to the worst 
one. 

Equation (Q) establishes a point-wise linear data model for our problem, that is, it holds true independently for each 
direction r or, likewise, for each pixel in the vector data map x. At this point, our problem can be formulated as follows: with 
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reference to the data model [Eq. (m)], for each pixel p identified by a unit vector f , estimate both the N- dimensional source 
vector s and the M x N mixing matrix A from the M-dimensional data vector x. 



3 THE FastICA ALGORITHM 

The problem of estimating both A and s from x via Eq. (Q) is unsolvable in the absence of additional information. The ICA 
approach assumes that 

• the components of vector s are independent random processes on the map domains; 

• the components of vector s, save at most one, have non-Gaussian distributions. 

In this section, we summarize an efficient strategy, described in detail in Hyvarinen & Oja (1997) and Hyvarinen (1999), to 
separate the independent components s, exploiting these assumptions. 

It has been proved (Comon 1994; Cardoso 1998) that if the components of s [Eq. (]?])] are independent, a solution to our 
problem can be obtained by searching a transformation W for the data vector x such that the transformed vector y = Wx 
has independent components. We should thus find a measure of independence among the components of y and maximize it 
for W. Intuitively, an equivalent approach to achieve independence is to maximize non-Gaussianity. Indeed, due to the central 
limit theorem, a variable resulting from a mixture of independent variables is "more Gaussian" than the original variables 
themselves. This means that finding a transformation such that the Gaussianity of the variables is reduced is equivalent to find 
a set of transformed variables that are "more independent" than the original ones. A strategy to find maximally non-Gaussian 
transformed variables is to select a suitable measure of non-gaussianity and then maximize it for the transform operator W. 
This approach can be shown to be equivalent to other theoretically sound approaches proposed in the literature (Comon 1994; 
Amari & Chichocki 1998; Cardoso 1998). Since Eq. (Q) is a noisy model, it is necessary to find a measure of non-Gaussianity, 
via a suitable defined neg-entropy, that allows us to adopt computing strategies that are robust against noise. 

An approximation to the neg-entropy function has been proposed by Hyvarinen & Oja (2000) and Hyvarinen (1999) as 
a measure of non-Gaussianity for the transformed variables and, if the noise has the features assumed in the previous section 
and its covariance matrix is known, the Gaussian moments of y are shown to be robust estimates of the desired function. It 
is possible to maximize non-Gaussianity by means of a Newton algorithm, after some preprocessing of the data. The rows of 
the transformation matrix W can be estimated one by one, and a noise-robust estimate is achievable. This does not mean 
that by linearly transforming x through W we obtain a robust estimate of s, especially in very noisy cases. To see this, with 
no loss of generality, we can assume M = N, and W as a robust estimate of A -1 . In this case, by applying W to our data, 
we obtain 

y = W(As + e) = s + A _1 e . (8) 

The noise term A _1 e can strongly affect the high-spatial-frequency components of this solution, and thus some regularization 
strategy must normally be adopted to separate the source signals once W has been estimated. We intentionally do not face 
this point in the present work, which is entirely dedicated to the study of the capabilities of FastICA by itself. 

Another general feature that one should keep in mind is that the quality of the estimated matrix W becomes worse 
when one or more signals are much stronger then the others, even if the reconstructed spatial patterns of the independent 
components are still good. Indeed, suppose for simplicity to have no noise, two frequency scalings and two independent 
components only. Then, from Eq. ([]) it can be easily seen that at each output yi one has 

y^ = (WiiAu + Wi2A 21 ) Sl + (WnA 12 + W a A 2 2)s 2 ■ (9) 

If, for example, si S> S2, then a good spatial reconstruction for which yi is a copy of si can be achieved even in the case where 
the coefficients of si and s 2 are comparable, or, in other words, the matrix W is not a robust estimator of A -1 . 

The algorithm needs a preprocessing step (Hyvarinen 1999), where the data maps are "quasi whitened". This step aims 
at reducing the number of unknowns with respect to the original problem. Let us suppose we know the covariance matrix 
£ of the system noise. The mean of the data at each frequency is removed (the offset of each independent component can 
be recovered in the end of the separation process, as we show below), and the covariance matrix C of the zero mean data is 
calculated by just evaluating the following expectation over the set of all available pixels: 

C = fi{xx T } . (10) 

The preprocessing procedure consists in evaluating a modified noise covariance and a quasi- whitened data set, respectively, 
as follows: 

£ = (C - £r 1/2 E(C - Er 1/2 , (11) 
x = (C - £)~ 1/2 x . (12) 
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The separation algorithm estimates the matrix W row by row, thus permitting to recover the components one by one. Let w 
be an M-dimensional vector whose dot-product by x gives a component of the transformed vector y = Wx (w T is a row of 
matrix W). In order to find a noise-robust estimation of one such vector, we apply the following iterative algorithm, equipped 
with a suitable convergence criterion, where g is a regular non-quadratic function whose form can be chosen as g(u) = u 3 , 
g(u) = tanh-u, g(u) = exp(— u 2 ) for most problems (see Hyvarinen 1999), and g is its first derivative: 

(i) Choose an initial vector w; 

(ii) update it through 

w nem = £{x<?(w T x)} — (I + t)wE{ g'(w T ±)} 

where again E denotes expectation over all the available samples; 

(iii) let 



(iv) Compare w new with the old one; if not converged, go back to (ii), if converged, begin another process. 

This procedure maximizes the non-Gaussianity of the component w T x. Let us assume that a certain number, k, of rows of 
W have been evaluated. To estimate the (k + l)-th row, we must search it in the subspace orthogonal to the first k rows. To 
achieve this goal, an orthogonalization step (e.g. Gram-Schmidt) must be inserted between steps (ii) and (iii). 

Once the separation matrix W has been found, physical normalization, frequency scalings and offsets of each independent 
component in the data can in principle be reconstructed. Consider normalization and frequency scalings first: by simply 
inverting Eq. (^|), we get 

x = W- 1 y, (13) 
which means that the independent component Xu present in x contributes to the signal at vt\i frequency: 
x\P (physical units at frequency u) = (W~ 1 )vjVj ■ (14) 
Then, the frequency scaling for the component j between frequencies v and v' is given by 

The offsets can be recovered as follows. True input data can be written as the sum of zero mean component plus mean values, 
indicated as x and (x) respectively. The zero mean component is reconstructed as in Eq. (|l3|), while the total signal obeys: 

(x)+x = W- 1 (y+(y» . (16) 

Then, by simply taking the mean of this expression we get 

(y> = W(x) , (17) 
so that, at i^th frequency, the offset of the independent component j is 

(x)l (physical units at frequency v) = {W~ 1 ) L .j{y)j . (18) 

Furthermore we can derive the signal-to-noise ratio in the reconstructed map. Since the system noise covariance matrix S is 
assumed to be known, it is possible to build noise constrained realizations n x for each frequency channel. Once the matrix 
W has been recovered, the corresponding noise realizations in the FastICA outputs n y are 

n y = Wn x . (19) 

The noise is then transformed like signals and an estimation of the noise into the reconstructed components is achievable. 
From these "processed" noise it is straightforward to obtain the signal-to-noise ratio. This is a useful tool when, e.g., the 
number of underlying components is not known a-priori. Suppose, for example, that the number M of observing channels is 
larger than the number N of the underlying sky components. Using all the channels leads to the recovery of M "components", 
the fictitious ones being characterized by a signal-to-noise ratio of the order, or below, unity. 



4 TOY SKY MODELS 

In the next Section we shall show FastICA results on simulated skies which are not too far from what the Planck satellite 
should observe, including the instrumental noise. However, we preferred to approach that case gradually, inserting here 
a Section which starts with highly idealized skies and ends just one step before the Planck case. We want to test the 
capabilities of the algorithm when components are well known and controlled, allowing a better understanding of the more 
realistic simulation of the next Section. 
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Table 1. Figures of merit for the reconstruction of the CMB, convolved with a 7° FWHM beam, and of a highly non-Gaussian signal 
(x 2 ) for different values of N s ^ e . The columns "moments" report the larger relative error on skewness and kurtosis. Simulations with 
(without) * have the x 2 signal 3 (30) times stronger than the CMB. Residual are evaluated for the CMB signal only. 
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Figure 1. FastICA reconstruction of the dust component (toy model II). The map is plotted in a non-linear scale to emphasize the 
high Galactic latitude structures, reflecting those of the CMB. 



4.1 Toy Model I 

In order to evaluate the ability of the FastICA algorithm to separate independent sky components, we start with a very 
simple toy model, comprising only two components (N = 2) with different distributions without, for the moment, taking into 
account the instrumental noise. This toy sky model is then "observed" at two frequencies, 53 and 90 GHz; therefore N = M. 
FastICA is a pointwise algorithm and we want to verify the dependence of the goodness of the reconstruction upon the total 
number of pixels in the maps, signal statistics and spatial distribution of the signal. We use sky maps pixelised according to 
the HEALPix scheme (Gorski et al. 1999) and we start working with N s id c = 64 corresponding to a pixel size of about 1°. Our 
sky emission model is composed of (i) a standard Cold Dark Matter CMB sky normalized to COBE (CDM, flat geometry, 
fif, = 0.05, 0.cdm = 0.95, h — 0.50, three massless neutrino families, scale-invariant Gaussian perturbations), which is a 
realization of a Gaussian process, smoothed with a 7° Gaussian beam, and (ii) a randomly distributed component (i.e. with 
no spatial features) with a x 2 distribution with one degree of freedom: 

x -l/2 e -x/2 

/x(x) = 2V2r(l/2) ' (20) 

F being the Gamma function. The latter is highly non-Gaussian with a kurtosis [Eq.(p3|) below] of £ 10000; no beam 

smoothing has been applied to it. For the CMB we adopt a black body spectrum at Tcmb = 2.726 K, while the x 2 signal is 
assumed to have a constant antenna temperature. We have run two cases in which the rms signal of the latter component is, 
respectively, ~ 30 and ~ 3 times larger than the CMB at 90 GHz. 

As mentioned in § [5] there are three kinds of non-linear functions that approximate the neg-entropy: we have verified that 
in most cases the function g(u) = u 3 works better than the other two and therefore we adopt it for the following analysis. 
The goodness of the reconstruction is quantified in terms of four figures of merit: 1) the correlation coefficient 
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Table 2. Same as in Table 1, but for a 1° FWHM beam. 
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Table 3. Figures of merit for the reconstruction of CMB and dust emission convolved with a Gaussian beam with 1° FWHM, for 
different Galactic latitude cuts. The t indicates a case in which the dust signal has been randomly distributed all over the sky. Residuals 
are evaluated only for the CMB map. 
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rs(x in ,x out ) = _ (21) 

yf E[(x in - X in ) 2 } y< E[{x out - X out ) 2 } 

between the output and input data x out ,Xi n (E again denotes the expectation value; x = E(x) represents the mean value), 
2) the percentage error in the frequency scaling of the reconstructed signal, 3) the percentage errors on higher order moments 

EUx-xf] , , 

Sk6WneSS = EKx-x) 2 )^ ' (22) 

E\<x-x) 4 } , . 

k-to S1 s=^-—ll-3 (23) 

of the reconstructed signal, and 4) the relative residual after re-calibration and subtraction of the corresponding input com- 
ponent. These figures of merit are evaluated also for a map with a higher number of pixels (A s id c = 128) and a higher angular 
resolution (smoothing with a 1° FWHM beam) of the CMB signal. As shown by Tables |l| and ^, the correlation between input 
and output signals is extremely tight for both components and does not vary as their intensity ratio changes by an order of 
magnitude. The correlation coefficient increases with increasing number of pixels. 

The recovery of the frequency dependence of each component is also good and obviously improves when the component 
becomes relatively stronger (while the recovery of the frequency spectrum of the other component is degraded roughly by the 
same factor). This result can be understood by consideration of Eq. (^): even when the spatial correlation is extremely good, 
the matrix W is not necessarily a robust estimator of A -1 . In the case of the weaker component, better accuracy is achieved 
for higher angular resolution and larger number of pixels. 

The very accurate recovery of higher order moments shows that the algorithm, at this level of simulation, does not 
introduce extra non-Gaussianity into the reconstructed signal. 

Equation (jf4|) allows us to convert the FastICA output into physical units. The residuals after subtraction of the re-scaled 
output map from the input one, show a pattern that is exactly the one of the other component: there is a sort of "mirroring" 
of one component into the other in the residuals. The level of residuals depends on iV s id e of the map; their rms values are 
about 10 -3 times smaller than the rms of the original map. 

For fixed A^ide, the quality of the reconstruction is significantly better in the case of the 1° resolution, thanks to the 
better sampling of the CMB statistics. 



4.2 Toy Model II 

We now replace the un-physical \ 2 signal with dust emission, modelled scaling to 53 and 90 GHz, with the spectral shape of 
Eq. (|27]), the maps by Schlegel et al. (1998). The CMB plus dust emission maps are then convolved with a Gaussian beam of 
1° FWHM. We adopt Aside = 128. The main difference with the previous model is that now we have inserted a signal (dust 
emission) with a strong spatial structure. We have considered this foreground, at frequencies where Galactic synchrotron and 
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Figure 2. Simulated skies at several Planck frequencies. Signals have been plotted in a non-linear scale to emphasize high Galactic 
latitude structures. 

free-free emissions should dominate, because its spatial distribution is known at high resolution (6') and we wanted to explore 
the effect of its markedly non-uniform sky distribution, highly concentrated on the Galactic plane. 

The power-law function, g(u) = u 3 , still works better than the others. A summary of the results, in terms of the usual 
figures of merit, is in Table |^. We note that, in the case of the all-sky analysis (first row of the Table), while the reconstructed 
CMB map is very well correlated with the input one, this is not the case for the dust. A significant degradation of the 
reconstructed dust map is clearly visible in Fig. [I], representing the output dust map, showing some CMB spatial features at 
high Galactic latitudes. Furthermore the derived frequency dependencies of both signal, but particularly that of the CMB, 
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Table 4. Skewness and kurtosis of dust emission and dust/CMB rms ratio for sky regions at low and high Galactic latitudes. 
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Table 5. Assumed parameters for the Planck channels considered here. Sensitivities, for a 14 months mission, are in antenna temperature. 
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deviate significantly from the expected one. In addition, when the CMB map is converted into physical units, a residual map 
showing dust spatial structures is present with an rms only 50 times smaller than the CMB rms, which is much worse than 
the result we obtained in Section 4.1. 

The main reason of these results is that the dust is highly not uniformly distributed across the sky. Indeed, as noted 
above, the recovery of the frequency dependence, i.e. of the separation matrix, is poor if one component is much weaker than 
the other; the situation worsens if the ratio of the two components is strongly varying across the map, as is the case here: dust 
dominates on the Galactic plane, while the CMB dominates at high Galactic latitudes. In addition, dust has different statistics 
on and off the Galactic plane, as Table. ^ shows. In order to check this point, we made an experiment separating CMB and an 
artificial "sky-uniform" dust map generated by randomly distributing its signal over all the sky (simply randomly shufnying 
the pixels so that the Galactic plane structure is destroyied): the quality of the reconstructed maps becomes extremely good 
(fourth row in Table |^) . 

These simulations show that a better reconstruction can be achieved for uniformly distributed sky signals. In the present 
problem, this can be obtained by analyzing separately regions around the Galactic plane and at high Galactic latitudes (second 
and third row of Table ^), where the Galactic emission has very different spatial distributions (see Table. |^). 



4.3 Effect of Instrumental Noise 

We now add instrumental noise to the toy model considered in the previous sub-section. The noise is assumed to be uniform 
across the sky and such that the average signal-to-noise ratio at high Galactic latitudes is ~ 2. Again, the all-sky reconstruction 
is unsatisfactory. However, if we analyze separately regions within |6| = 25° and at higher Galactic latitudes, we recover the 
CMB frequency scaling with a 0.5% error at high Galactic latitudes and that of dust emission with the same error at low 
latitudes. The signal-to-noise ratio at high Galactic latitudes is ~ 2 for the CMB (as in the input map) and ~ 1 for dust. At 
low Galactic latitudes the signal-to-noise ratio is still ~ 2 for the CMB, and ~ 12 for dust emission. 



5 FastICA AND PLANCK SIMULATED MAPS 



We simulate all-sky maps adopting the HEALPix pixelization scheme (Gorski et al. 199£ ) . The maps are in antenna tempera- 
ture and, as described in detail below, include thermal dust, free-free, and synchrotron emissions, on top of CMB fluctuations. 
For the CMB we assume a black body spectrum with Tcmb = 2.726 K, so that fluctuations in antenna temperature, STa,cmb 
are related to brightness temperature fluctuations STcmb by: 



x 2 e- 



<52a,cmb(^) = 7— tt^STcmb , (24) 

[e x — iy 

where x = hv/kTcMB- We have adopted brightness temperature fluctuations, STcmb, corresponding to a flat ACDM cos- 
mology with Qcdm = 0.35, Qb — 0.05, £7a = 0.6 and h = 0.65, three massless neutrino families and an initial scale-invariant 
Gaussian spectrum. Galactic foreground emissions are described in detail below. 

The CMB and foreground maps are then summed up together with simulated white noise maps, generated by assuming 
the rms noise fluctuation level expected for the Planck channels (see Table |^). The simulated skies can be seen in Fig. [| 
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Figure 3. Input (left) and output (right) synchrotron maps for case 1, at 30 GHz. Maps have been plotted in a non-linear scale to 
emphasize high Galactic latitude features. 



Given our sky composition, three cases exploiting three subsets of Planck frequency channels are considered. The first 
case (case 1) makes use of the two lowest frequency channels (30 and 44 GHz) assuming a mixture of CMB plus the dominant 
diffuse Galactic foreground, which is assumed to be the synchrotron emission. Note that a significant contribution from 
Galactic free-free should be present in the real sky, at least at low Galactic latitudes, but a good template for it is still 
missing. For the moment, we neglect it. We apply FastICA to all-sky maps as well as to regions of high (|6| > 20°) and low 
(|6| < 20°) Galactic latitude. 
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Table 6. Percentage errors on FastICA frequency scaling reconstruction. 



Considered frequencies (GHz): 


30/44 


100/143 


143/217 


217/353 


Case 1, all-sky 


CMB 

Synchrotron 


1 

3 








Case 1, |6| > 20° 


CMB 

Synchrotron 


1 

0.3 








Case 1, |6| < 20° 


CMB 

Synchrotron 


2 
3 








Case 2, all-sky 


CMB 

Thermal component 
Synchrotron 




1 

0.4 
75 


0.6 
0.2 

84 




Case 3, all-sky 


CMB 

Thermal component 








46 
0.3 


Case 3, |6| > 20° 


CMB 

Thermal component 








12 
0.7 


Case 3, |6| < 20° 


CMB 

Thermal component 








64 
0.03 



The second application (case 2) is more CMB oriented: it refers to the Planck "cosmological" channels at 100, 143 and 
217 GHz, where foregrounds are at a level comparable with noise; in this case, we apply FastICA only to all sky maps. 

The third case (case 3) considers again two signals and two channels, but we move to higher frequencies, 217 and 353 
GHz, where the CMB and our thermal Galactic component, almost coincident with dust, dominate over synchrotron, which 
was neglected. Again we apply FastICA to all sky maps, to the region \b\ < 20° and to the complementary region at high 
Galactic latitudes. 

Before to go to the results, let us describe briefly the methodology used. Simulated maps were originally produced with 
pixels of size 3.5', corresponding to N s id e = 1024 in the HEALPix pixelization scheme, and such high- resolution maps should 
be taken as "reference" skies. However, we can note from Table ^| that the highest frequency channels have 5' resolution, which 
in principle requires a smaller pixel size to Nyquist sample the beam-width. Since in this preliminary study we are neglecting 
the astrophysical signals relevant on the smallest scales (such as extragalactic point sources), we restrict our analysis to 
multipoles £ ,$ 2000, for which pixels of 3.5' are enough. 

The reference maps have been smoothed with the appropriate FWHM of Planck channels, and the corresponding 
instrumental noise was added. Due to Assumption 3 in Sect. B the FastICA algorithm requires input maps having the same 
angular resolution. This situation has been achieved by reducing, for each case considered, the effective working angular 
resolution FWHM wor k to the worst one. This implies an extra-smoothing of the observed (signal + noise) maps having 
FWHM<F WHM W ork. Note that the Planck nominal noise (where the noise is assumed uniformly distributed on the sky) is 
generally given per resolution element corresponding to beams of different size, depending on the channel. This extra-smoothing 
process required by FastICA changes the noise rms with respect to the value reported in Table |E| the new rms values are 
used to construct the £ matrix in Eq. ([H]). Simulations on cases 2 and 3 have been performed maintaining N ai d e = 1024, 
while in the case 1 we worked with N S ide = 256 because of the smaller angular resolution of the lowest Planck channels. 

We quantify the quality of FastICA outputs with reference to the percentage errors on frequency scaling recovery, 
on normalization, and on the angular power spectrum. For each application we estimated the relative error on the scaling 
reconstruction between nearby frequencies, by comparing the ratios of the W _1 matrix elements in Eq. (^) with their input 
values. The results, for all the cases, are shown in Table []. The normalization allows us to measure in some detail how the 
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Figure 4. Input (left) and output (right) CMB maps for case 1, at 30 GHz. 



noise affects the separation matrix, since as we have seen in § 2, Eq. (|l4[), the normalization of the FastICA outputs to 
physical units involves a direct product of an element of the inverted separation matrix by a FastICA output. The angular 
power spectrum of each physical component is estimated using the FastICA output maps as follows. At each frequency, and 
for each signal, the observed Ce are the sum of signal and noise contributions: 

C™ ap = C e W{£) + C"°" e 

where the filter function 
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W(£) = exp[-£(£ + l)(FWHM/rad) 2 /(21og2)] (25) 

describes the effect of the beam shape on multipoles and C" ms,! represents the noise contamination, which is additive only 
if noise and signal are perfectly uncorrelated. Our guesses about the noise rms allow us to Monte Carlo simulate the C" mse . 
The power spectra of FastICA output noisy maps, obtained with the code Anaf ast (HEALPix package), should therefore be 
a good estimate of the quantities [CtW(£) + C"° lse ] for the signals observed. We can invert the previous relation and estimate 
the power spectra of the signals, since we know W(£) and have a guess for C™° !se : 



Cmap /-inoi, 
e ~ ^ P 



(26) 



W(£) 

This is the relation used in the following analysis, to evaluate the performances of FastICA separation in terms of the angular 
power spectrum. 



5.1 Case 1 

Components to be recovered are CMB and synchrotron which are assumed to dominate the sky emission at 30 and 44 GHz. 
The synchrotron emission template is the 408 MHz map of Haslam et al. (1982), extrapolated to the considered frequencies 
assuming a uniform spectral index (3 syn — —2.9. On account of the poor resolution of the Haslam map (~ 0.85°) we artificially 
add Gaussian small scales fluctuations with a power spectrum Pi tx £~ 3 . 

As already mentioned, we have applied the FastICA algorithm to the whole sky, as well as to regions around the Galactic 
plane (|6| < 20°) and of high Galactic latitude (|6| > 20°). The 44 GHz noisy map was smoothed to a FWHM work = 33.6', 
corresponding to the Planck resolution at 30 GHz. We finally re-gridded the maps with pixels of size ~ 14', (corresponding 
to iV s id c =256), which is enough to Nyquist sample the working beam size. Figures show input versus output maps, while 
figures ^,^| give the corresponding power spectra, showing, from top to bottom, the all-sky, high and low Galactic latitude 
results, respectively. 

Synchrotron is reconstructed remarkably well over whole sky. The residual noise shows up in the high £ tail of the 
reconstructed spectrum. The normalization is also correctly recovered. The reconstruction for regions at low and high Galactic 
latitudes is still good, although it does not reach the quality of the all-sky analysis. Note however that the synchrotron map has 
essentially no power above £ ~ 100; this is the reason why we plotted this signal only up to multipoles £ ~ 300 corresponding 
to scales slightly below 1°. When the angular power spectrum of this signal is estimated for a fraction of the sky, oscillations 
are found, which propagate down to very low £s. Such oscillations have been smoothed by averaging over intervals A£ — 10. 
However FastICA has recovered the correct normalization of the reconstructed synchrotron, as can be seen by looking at the 
scales in Figs. ^| and ^. 

The CMB is reconstructed well on all the scales down to the adopted resolution scale, which corresponds to £ ~ 500. 
Also the normalization was recovered quite correctly. The rise in the recovered spectra at £ ~ 500 is because the maps have 
been smoothed and then gridded before separation, as we explained above: this introduces extra power on the spectrum, at 
multipoles corresponding to the pixel size, which is enhanced by dividing by the beam window function as in Eq.(j2(j). 

The percentage errors on the reconstructed frequency scalings, reported in Table ^j, are at the percent level for both 
components. 



5.2 Case 2 

In this frequency range (100-217 GHz) the CMB is the dominant component, except for the Galactic plane region. Our 
simulated skies contain a mixture of CMB, synchrotron and thermal emission, described below. The dust emission templates 
have been obtained by extrapolating the maps by Schlegel et al. (1998), which combine IRAS and DIRBE data, assuming a 
grey-body spectrum (expressed in antenna temperature) , 

i> p+1 hv 

(27) 

with uniform temperature T^ust = 18 K and emissivity f3 — 2. 

Unfortunately, no maps of free- free emission are available at the moment, although an all-sky map of H a emission (which 
is known to be a good tracer of free-free emission) will be available in a couple of years (Reynolds & Haffner, 2000). On 
the other hand, at ^ 100 GHz free-free emission may exceed the dust emission component. In order to simulate the free-free 
contribution we assume, somewhat arbitrarily, that it is perfectly correlated with the dust itself, i.e. that it has the same 
spatial distribution. Its antenna temperature scales with frequency as Taj/ oc v~^ff , with /3// = —2.1. The relative amplitude 
of dust and free- free emission is assumed to be a factor of 3 at 100 GHz (De Zotti et al. 1999). Thus, for the FastICA analysis, 
free-free and dust emission effectively behave as a single component, referred to, in the following, as the "thermal" component, 
having a spectrum described by: 
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synchrotron power spectrum at 30 GHz, all — sky aralysis 



FastICA all — sky reconstructia n 




synchrotron power spectrum at 30 GHz, off Galactic plane analys 



FastICA off Galactic plane reconstruction 



300 



Figure 5. Input (left) and output (right) synchrotron angular power spectra for case 1, at 30 GHz. From top to bottom we plotted 
all-sky, low and high Galactic latitude FastICA results, respectively. 



fthermaliv) - 3 ( 100GHz ) 



V \' 3ff , fdust(v) 



/ dus t(100GHz) • (28) 

The Galactic signals are well below the noise at medium and high Galactic latitudes so that no reconstruction is possible 
there at these frequencies. Therefore we have restricted the analysis to all-sky maps. The smoothing of the noisy sky maps 
has been performed with an FWHM of 10', and we have used iVside=1024 (3.5' pixels). Figures (?) and ^ show the input (left) 
and output (right) maps and power spectra, respectively, at 100 GHz. 

The CMB power spectrum is reconstructed up to i ~ 1500 with the correct normalization. On the contrary, thermal 
emission and synchrotron power spectra are much more affected by noise. The power spectrum of thermal emission is quite 
well reconstructed up to I ~ 600, but its normalization is overestimated by a factor ~ 1.8. The shape of the synchrotron 
power spectrum is recovered up to I ~ 100, but its normalization was lost. 
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Figure 6. Input (left) and output (right) CMB angular power spectra for case 1, at 30 GHz. From top to bottom: all-sky, low, and high 
Galactic latitude FastICA results, respectively. 



The reason of the bad quality of the synchrotron reconstruction resides in the weakness of this component at these 
frequencies, with respect to noise, and in the poor statistics contained in the synchrotron template we use, as we already 
stressed. 

However, it is remarkable that FastICA could separate, to some extent, thermal and synchrotron emissions, despite their 
low signal-to-noise ratio and the correlation in their spatial distributions (both are maximum on the Galactic plane). The 
latter point illustrates the fact that the statistical independence required for FastICA to work is not so much related to the 
spatial pattern of the components to be separated but, rather, to their distribution functions. The templates we have adopted 
for thermal and synchrotron emissions do have quite different distribution functions, as illustrated by Table [j], giving their 
skewness and kurtosis. This makes FastICA able to perform some separation of them. 
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Figure 7. Input (left) and output (right) maps for case 2, at 100 GHz. 



As we can see from Table g frequency scalings are recovered very well, at the percent level, for the CMB and thermal 
components, while the error is considerably larger for synchrotron, which is the weakest component in this case. 

5.3 Case 3 

For this analysis, we selected two channels, 217 and 353 GHz, and two components, CMB and thermal emission which, at 
these frequencies, is dominated by dust. Maps with 7V s id o =1024 have been smoothed with a FWHM of 5'. As in case 1, we 
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Figure 8. Input (left) and output (right) angular power spectra of the CMB, of the thermal component, and of synchotron emission 
case 2, at 100 GHz. 



Table 7. Skewness and kurtosis of CMB, thermal component and synchrotron templates. 



Component 


Skewness 


Kurtosis 


CMB 


0.014 


0.087 


Thermal 


51.7692 


5098.16 


Synchrotron 


7.717 


80.693 
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Figure 9. Input (left) and output (right) CMB maps for case 3. 



have made separate analyses for the whole sky, and for the high and low Galactic latitude regions. The average signal-to-noise 
ratio for the dust is 190 over the whole sky, 246 for \b\ < 20°, and 6.6 at high Galactic latitude; for the CMB it is ~ 5.7. 

The results are reported in Figures, [u], [HJ and referring to 217 GHz. The shapes of the power spectra of both 
signals are reconstructed very well up to £ ~ 2000 in all cases shown. The amplitude of the power spectrum of thermal emission 
is accurately recovered in all cases, while in the case of the CMB this happens only at high Galactic latitudes, where the two 
components are comparable. As shown by Table the frequency scaling of thermal emission is recovered to a fraction of 1%, 
while for the CMB the error is much larger. 
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Figure 10. Input (left) and output (right) thermal emission maps for case 3, at 217 GHz. 



6 DISCUSSION AND CONCLUSION 

We have presented a first application of a new, computationally fast, technique for separation of astrophysical components, 
based on Independent Component Analysis concepts, FastICA. The technique has been applied both to toy models, useful 
to highlight some of its key properties, and to simulated all-sky maps, with reference to the forthcoming Planck mission 
(Mandolesi et al. 1998; Puget et al. 1998). 

We have achieved a considerable improvement over the previous application of ICA-type algorithms to astrophysical 
component separation (see Baccigalupi et al. 2000) . First of all, the method is considerably faster: processing of all-sky maps 
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Figure 11. Input (left) and output (right) CMB angular power spectra for case 3, at 217 GHz. From top to bottom: all-sky, low and 
high Galactic latitude results, respectively. 



with ~ 10 7 sky pixels (pixel size of 3'.5) took approximately ten minutes on a Pentium III 600 MHz computer. This is much 
less than required by any existing non-blind separation methods like Maximum Entropy, taking 6 hours on 30 processors at 
equivalent resolution (Stolyarov et al. 2001). Second, we have shown that FastICA works in the presence of instrumental noise, 
although under the simplifying assumption that it is white and uniformly distributed. Third, the algorithm was successfully 
applied both to all-sky maps and to limited portions of the sky. Using the full all sky maps obviously minimizes the sample 
variance. On the other hand, we have shown that the foreground removal improves if we cut regions (around the Galactic 
plane) where their distribution is highly non-uniform. 

Furthermore, we extended the formalism of the separation algorithm in order to obtain the normalization constants 
to convert FastICA output maps in physical units. We found that the accuracy on this calibration is limited only by the 
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Figure 12. Input (left) and output (right) angular power spectra of the thermal component for case 3 at 217 GHz. From top to bottom: 
all-sky, low, and high Galactic latitude results, respectively. 



instrumental noise. The same noise effect limits the frequency scaling reconstruction. On the other hand, FastICA still 
requires that all the input maps have the same angular resolution. 

We applied FastICA to different combinations of Planck frequency channels, namely 30, 44, 100, 143, 217, 353 GHz, 
containing a mixture of CMB, synchrotron and thermal emissions, the latter comprising free-free plus thermal dust emission. 
The nominal average instrumental noise levels and angular resolutions of Planck have been adopted. 

With all these ingredients, FastICA is able to obtain an excellent reconstruction (well calibrated) of all these components. 
The angular power spectra of both CMB and thermal emissions can be accurately reconstructed up to multipoles £ ~ 2000. 
Synchrotron is also recovered on all scales where our template has significant power, i.e. up to I ~ 300. Frequency scalings 
and normalization are recovered up to better than 1% precision. 
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Nevertheless, a lot of work to assess the real capabilities of this and other component separation algorithms for CMB 
experiments has still to be done. Here we list our main simplifying assumptions, which we have to relax in future. 

The spectral properties of synchrotron and dust emissions have been assumed to be constant over the sky, although it 
is known that both the synchrotron spectral index and the dust temperature vary from place to place. The spatial pattern 
of free-free emission was assumed to follow that of dust, which is certainly over-simplistic. Spinning dust emission was not 
taken into account. We also neglected extra-galactic sources, which are expected to contribute significantly on all Planck 
channels, as well as the thermal and kinetic Sunyaev-Zeldovich effect. Very promising point-source extraction techniques are 
being extensively studied (Cayon et al. 2000; Vielva et al. 2001). 

Finally, the instrumental noise is not uniform across the sky, as assumed here, but has a complex pattern determined by 
the adopted scanning strategy adopted. Also the instrumental beam is not perfectly circularly symmetric, etc. 

Despite all these still unaddressed issues, we demonstrated that FastICA is a very promising technique because of three 
main features: (i) the good separation achieved on our simulated maps, which contained noise at the Planck nominal level, 
(ii) the high computational speed, and (Hi) the extreme flexibility of the code, which is able to work on all-sky maps as well 
as on arbitrarily shaped sky regions in a straightforward way. We recall also that FastICA is a blind separation algorithm, 
able to estimate not only independent components present in the data, but also their frequency scalings. 
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